data{
int<lower=1> N; //number of observations
int<lower=1> C; //number of cities
int<lower=1> K; //nubmer of covariates
int<lower=1> P; //number of parites
int<lower=1, upper=C> city[N];
int<lower=1, upper=P> party[N];
matrix[N,K] x;
vector<lower=0, upper=1>[N] y;
}

transformed data{
matrix[N,K] Q_ast=qr_Q(x)[,1:K]*sqrt(N-1);
matrix[K,K] R_ast=qr_R(x)[1:K,]/sqrt(N-1);
matrix[K,K] R_ast_inverse=inverse(R_ast);

}


parameters{
vector[K] beta_tilde;
vector[C] ranef_c;
vector[P] ranef_p;
real<lower=0, upper=100> sigma_p;
real<lower=0, upper=100> sigma_c;
real<lower=0, upper=100> sigma_y;
real alpha;
real<lower=1, upper=7> nu;

}

transformed parameters {
vector[K] beta=R_ast_inverse*beta_tilde;
vector[N] y_hat;
vector[N] eta=Q_ast*beta_tilde;

for(n in 1:N){
y_hat[n]=alpha+eta[n]  + ranef_p[party[n]] + ranef_c[city[n]] ;
}


}

model{
//priors
target+=normal_lpdf(alpha|0,1);
target+=student_t_lpdf(sigma_c|4,0,1);
target+=student_t_lpdf(sigma_p|4,0,1);
sigma_y~exponential(1);
target+=gamma_lpdf(nu|2,0.1);
//likelihood


for(c in 1:C){
  target+=normal_lpdf(ranef_c[c]|0, sigma_c);
}

for(p in 1:P){
  target+=normal_lpdf(ranef_p[p]|0, sigma_p);
}



for(n in 1:N)
target+=student_t_lpdf(y[n]|nu, y_hat[n], sigma_y);

}

//for(n in 1:N)
// target+=normal_lpdf(y[n]| y_hat[n], sigma_y);

//}

generated quantities{
vector[N] y_rep;
vector[N] log_lik;
vector[K] fd_beta;



for(n in 1:N){
log_lik[n]=student_t_lpdf(y[n]|nu, y_hat[n], sigma_y)  ;
}

for(n in 1:N){
y_rep[n]=student_t_rng(nu, y_hat[n], sigma_y);
}

//for(n in 1:N){
//log_lik[n]=normal_lpdf(y[n]|y_hat[n], sigma_y)  ;
//}

for(n in 1:N){
y_rep[n]=normal_rng(y_hat[n], sigma_y);
}


for(k in 1:K){
  if(k==4)
  fd_beta[k]=(alpha+beta[k]*1)-(alpha+beta[k]*0);
  else if(k==5)
  fd_beta[k]=(alpha+beta[k]*1)-(alpha+beta[k]*0);
  else
  fd_beta[k]=(alpha+beta[k]*(mean(x[,k]+sd(x[,k]))))-(alpha+beta[k]*(mean(x[,k]-sd(x[,k]))));
}



}
